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CN ■ Combining conservation of energy throughout nearly-spherical collapse of 

Ph 1 galaxy clusters with the virial theorem, we derive the mass-temperature rela- 

C/} ■ tion for X-ray clusters of galaxies T = CM 2 / 3 . The normalization factor C and 

00 1 the scatter of the relation are determined from first principles with the additional 

assumption of initial Gaussian random field. We are also able to reproduce the 
CN ■ recently observed break in the M-T relation at T ~ 3 keV, based on the scatter 

q ■ in the underlying density field for a low density ACDM cosmology. Finally, by 

combining observational data of high redshift clusters with our theoretical for- 

O ■ 

■ malism, we find a semi-empirical temperature-mass relation which is expected to 

hold at redshifts up to unity with less than 20% error. 
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1. Introduction 

The abundance of clusters of galaxies provides one of the strongest constraints on cosmo- 
logical models (Peebles, Daly, & Juszkiewicz 1989; Bahcall & Cen 1992; White, Efstathiou, 
& Frenk 1993; Viana & Liddle 1996; Eke, Cole, & Frenk 1996; Oukbir, Bartlett, & Blan- 
chard 1997; Bahcall, Fan, & Cen 1997; Pen 1998; Cen 1998; Henry 2000; Wu 2000) with an 
uncertainty on the amplitude of density fluctuations of about 10% on ~ 10/i _1 Mpc scale. 
Theoretically it is often desirable to translate the mass of a cluster, that is predicted by 
either analytic theories such as Press-Schechter (1974) theory or N-body simulations, to the 
temperature of the cluster, which is directly observed. Simple arguments based on virial- 
ization density suggest that T oc M 2//3 , where T is the temperature of a cluster within a 
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certain radius (e.g., the virial radius) and M is the mass within the same radius. However, 
the proportionality coefficient has not been self-consistently determined from first principles, 
although numerical simulations have frequently been used to calibrate the relation (e.g., 
Evrard, Metzler, & Navarro 1996, hereafter EMN; Bryan & Norman 1998; Thomas et al. 
2001). 

It is noted that the results from different observational methods of mass measurements 
are not consistent with one another and with the simulation results (e.g., Horner, Mushotzky, 
& Scharf 1999, hereafter HMS; Neumann, & Arnaud 1999; Nevalainen, Markevitch, & For- 
man 2000, Finoguenov, Reiprich, & Bohringer 2001, hereafter FRB ). In general, X-ray mass 
estimates are about 80% lower than the predictions of hydro-simulations. Fig 1 compares 
X-ray cluster observational data with best fit line to EMN simulation results (FRB). On 
the other hand, mass estimates from galaxy velocity dispersion seem to be consistent with 
simulation results (HMS). The error in the gravitational lensing mass measurements is still 
too big to distinguish between these two (Hjorth, Oukbir & van Kampen 1998). 
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Fig. 1. — Observational mass-temperature data based on temperatures of X-ray clusters 
(crosses with la error bars) vs. the best fit line to EMN simulation results (FRB). 
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Another recent observational finding is the possible existence of a break in the T — M 
relation. By use of resolved temperature profile of X-ray clusters observed by ASCA, FRB 
have investigated T — M relation in the low-mass end and find that M oc T~ 2 , compared 
to M oc T~ 3 / 2 at the high mass end. Suggestions have been made to explain this behavior 
by attributing it the effect of formation redshift (FRB), cooling (Muanwong et al. 2001) and 
heating (Bialek, Evrard & Mohr 2000) processes. 

In this paper we use conservation of energy for an almost spherically collapsing region 
to derive the M-T relation. In §2 we find the initial and final energy of the cluster. §3 
constrains various factors which enter the normalization of M-T relation, via statistical 
methods, simulation and observational input. §4 considers predictions of our model and 
comparison with observational and simulation results. In §5, we discuss the limitations of 
our approach and justify some of the approximations. §6 concludes the paper. 



2. Conservation of Energy 

2.1. Initial Kinetic and Potential Energy of a Proto-Cluster 

We begin by deriving the kinetic energy of the proto-cluster. We can write velocity as 
a function of gravitational potential (pi (e.g. Padmanabhan 1993): 

v = ^x--|-V^, (1) 

at the initial time in the linear regime, where Hi is the Hubble constant at the initial time. 
There is a small dependence on the initial density parameter fli in equation (1) which we 
ignore since initially it is very close to unity and the difference will be in second order terms 
that we ignore for the proto-cluster. Then the kinetic energy is given by: 

K i = \J P v2d " x = /(! + 5 i)\ H ^ ~ ^V0,| Vx. (2) 
Keeping the terms up to the linear order we obtain 

Ki = X - Pl J (Hfx 2 + ^V%x 2 - ^x.V0,)rf 3 x. (3) 
In deriving equation (3), we have used the Poisson equation: 

V 2 0, = 47rGM, (4) 



to substitute for 5 iy where pi is the initial mean density of the universe. For equation (3) we 
then use Gauss theorem to make the third term similar to the second, at the expense of a 
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surface term: 



K t = X - Pl J (Hfx 2 + ^x 2 V 2 fa)d 3 x - X - Pi j> x 2 Vfa.d*. (5) 

Assuming that the deviations from spherical symmetry is not important at the boundary of 
the proto-cluster, we find Vfa in the second term in equation (5) as a function of Sf. 



„ GSM ^G Pi f . l3 



where G is the gravitational constant and Ri is the boundary radius of the initial proto- 
cluster, which leads to 

Kt = J [(1 + 25 l )x 2 - R%\d\ (7) 

where we have used equation (4) to substitute back for fa and also the definition of f2j = 
8nG Pl /(3H 2 ). 

Let us now find an expression for the gravitational potential energy of the proto-cluster. 
Using its definition we have 

Uj = _M[ A(i^,fa))(i^.W) |# , iA2 . (8) 

2 J J \Xi — X2\ 

Keeping the terms to the first order and using the symmetry under interchange of x\ and 
x 2 , we arrive at 

G Pi , or /_ \nj3_ f d ^ x 2 



Ui 



2 

which, taking the second integral in a spherical volume, gives 



f(l + 2S i (x 1 ))d 3 x 1 [ dX2 (9) 
J J \x 1 -x 2 \ 



Ui = ~fii [(1 + 25 l ){^-^)dFx. (10) 



Adding equations 7 and 10 gives the total initial energy 3 

4nG Ait 



Ei = 



3 

to the first order. Defining x, Si and B as: 



(1 - Sli)R\ - \ J 6 % {x){R 2 - x 2 )d*x], (11) 



5 = -^-, 5 i = 5 i + |(f2 i _l),B= / 6i{x)(l-i?)<Px, (12) 



3 It is easy to see that 1 — = 0(5») and so we have neglected it in the first order terms. 
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equation (11) is simplified: 

E, = -^f^B. (13) 

The integral in the definition of B is in fact a three dimensional integral and the limits denote 
that the integration domain is the unit sphere. Note that Si can be considered as the density- 
perturbation to a flat (i.e. Q tot = 1) universe for which energy vanishes. Another aspect 
of this statement is that both terms in the definition of 5^ scale as a in the early (matter 
dominated) universe and so does 5i itself. For a flat universe the first term dominates at 
high redshift since the second term scales as a 3 . 



2.2. Energy of a Virialized Cluster 

According to the virial theorem, the sum of the total energy Ef of a virialized cluster 
and its kinetic energy Kf vanishes. However, non-vanishing pressure at the boundary of the 
cluster can significantly modify the virial relation. Integrating the equation of hydrostatic 
equilibrium, we have 

K f + E f = 3P ext V, (14) 

where P ext is the pressure on the outer boundary of virialized region (i.e. virial radius) and 
V is the volume. For now we assume that the surface term is related to the final potential 
energy Uf by 

3P ex tV = -uU f . (15) 

We will consider the coefficient v and its possible mass dependence in §3.4. For a system of 
fully ionized gas plus dark matter, the virial relation (14) with equation (15) leads to 

+ rr 3,, 9 3M aas kT . . 

- i — )E/ = K / = -M DM vt + ^-, (16) 

where a v is the mass-weighted mean one-dimensional velocity dispersion of dark matter 
particles, M^m is the total dark matter mass, k is the Boltzmann constant, fi = 0.59 is the 
mean molecular weight and m p is the proton mass. Assuming that the ratio of gas to dark 
matter mass in the cluster is the same as that of the universe as a whole and / is the fraction 
of the baryonic matter in the hot gas, we get 

(17) 

with P S pec = cr^ / (kT / fxrrip) . Hydrodynamic simulations indicate that f3 spec ~ 1. For simplic- 
ity we define f3 spec as 

/Vc = /Wl + (fP7pec - 1)^]- (18) 
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So equation (16) reduces to: 



3P spec MkT 

L f 



Assuming energy conservation (i.e., Ei = Ef) and combining this result with equations (13, 
16) lead to the temperature as a function of initial density distribution: 

kT = ^ ( i±ji )HMB (20) 

In the next subsection §2.3 we will find an expression for HfRf in terms of cluster mass M 
and the initial density fluctuation spectrum. 



2.3. Virialization Time 

Defining e as the energy of a test particle with unit mass, which is at the boundary of 
the cluster Ri of mass M initially: 

e=*-™. (21) 
2 R, 1 ' 

The collapse time t can be written as 

2irGM 



(-2e) 



3 ■ 

2 



(22) 



Following the top-hat model, we assume that the collapse time of the particle is approxi- 
mately the same as the time necessary for the particle to be virialized. Choosing t to be the 
time of observation and assuming the mass M, interior to the test particle is virialized at t, 
and by combing equations (1, 6, 22), we find e as function of initial density distribution and 
relate it to the collapse time: 

-2e = AflfiJ? jf Ux)d*x = (^)f . (23) 

Using M = ^irR^pi and the Friedmann equations we obtain 

^jf*^ = !<&'■ (24) 
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2.4. M-T Relation 

Combining equations (20) and (13,24), we arrive at the cluster temperature-mass rela- 
tion: 

4 T = (i^)(i±^)«§(f). (25) 

Notice that although B and A are both functions of the initial moment t i: since both are 
proportional to the scale factor a, the ratio is a constant; the derived T-M relation (equation 
25) does not depend on the adopted initial time, as expected. As a specific example, in the 
spherical top-hat model in which the density contrast is assumed to be constant, this ratio 
B/A is |. Let us gather all the unknown dimensionless factors in Q: 

Q = (^)-\\^)(j)(Hty 2/3 - (26) 

Then, inserting in the numerical values, equation (25) reduces to: 

kT = (6.62 k eV)Q (l5I ^P, (27) 

or equivalently: 

kT 

M = 5.88 x 10 13 Q- L5 (— — ) L5 /i- 1 M (28) 

i Ke v 

where H = lOO/i km/s/Mpc is the Hubble constant. This result can be compared with the 
EMN simulation results: 

kT 

M 200 = (4.42 ± 0.56) x lO^^^) 1 - 5 ^^. (29) 

To convert the M 50 o masses of EMN to M 2 oo we have used the observed scaling of the 
mass with density contrast Ms oc 5~ ' 266 (HMS), which is consistent with the NFW profile 
(Navarro, Frenk, & White 1997) for simulated dark matter halos as well as observations 
(e.g., Tyson, Kochanski, & dell' Antonio 1998), in the relevant range of radius. 



3. Numerical Factors, Scatter and Uncertainties 

In this section we try to use different analytical methods as well as the results of dark 
matter simulations and observed gas properties, available in the literature, to constrain the 
numerical factors which appear in Q (equation 26). 



- 9- 



3.1. (3 spec and the gas fraction 

So far we have made no particular assumption about the gas dynamics or its history, 
4 and so we are going to rely on the available observational results to constrain gas prop- 
erties. (3 sp ec is defined as the ratio of kinetic energy per unit mass of dark matter to the 
thermal energy of gas particles. This ratio is typically of the order of unity, though different 
observational and theoretical methods lead to different values. 

The hydrodynamic simulation results usually point to a larger value of f3 spec . For ex- 
ample, Thomas et al. (2001) find (3 spec = 0.94 ± 0.03. On the other hand, observational 
data point to a slightly lower value of I3 spec . Observationally there is yet no direct way of 
accurately measuring the velocity dispersion of dark matter particles in the cluster and one 
is required to assume that the velocity distribution of galaxies follows that of dark matter or 
adopt a velocity bias. Under the assumption of no velocity bias Girardi et al. 1998 find it to 
be 0.88 ± 0.04. Girardi et al. 2000 study (3 S p ec for a sample of high redshift clusters and do 
not find any evidence for redshift dependence. From the theoretical point of view, the actual 
value of f3 spec might be substantially different from the observed number, because both the 
velocity and density of galaxies do not necessarily follow those of the dark matter, which 
could have resulted in some non-negligible selection effects. Unknown sources of heating such 
as due to gravitational energy on small scales which is often substantially underestimated 
in simulations due to limited resolutions or baryonic processes like supernova feedback may 
affect the value of f3 spec as well. 

Hydrodynamic simulations show that only a small fraction of the baryons contribute 
to galaxy formation in large clusters (e.g., Blanton et al. 2000) and so / is close to one. 
Observationally we quote Bryan (2000) who compiled different observations for the cluster 
mass fraction in gas and galactic component with 

/ = 1 - 0.26(T/10 keV)~ a35 , (30) 

albeit with a large scatter in the relation. Inserting equation (30) into equation (18), we see 
that the correction to j3 spec ~ 0.9 is less than 5% for all of feasible cosmological models which 
are dominated by non-baryonic dark matter. In what follows, unless mentioned otherwise, 
we adopt the value f3 spec = 0.9 and absorb any correction into the overall normalization of 
the T-M relation. 



4 With the exception of any heating/cooling of the gas being negligible with respect to the gravitational 
energy of the cluster. 
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3.2. B/A: Single Central Peak Approximation and Freeze-Out time 

In the original top-hat approximation (Gunn & Gott 1972) which has been extensively 
used in the literature, the initial density distribution is assumed to be constant inside Ri, 
which leads to B/A = 2/5. Shapiro, Iliev, & Raga (1999) and Iliev & Shapiro (2001) have 
extended the original treatment of top-hat spherical perturbation to a more self-consistent 
case with a trancated isothermal sphere final density distribution including the surface pres- 
sure term (see §3.4). Here we consider the general case with an arbitrary density profile of 
a single density peak. But before going further, let us separate out the term due to space 
curvature in the definition of B/A. Using the definitions of B and A in equations (12) 
and (24), and Friedmann equations to insert for pi in terms of the present day cosmological 
parameters, we get: 

B b ,b 2 W n _ w Ht x 2 

A a a 5 J£ m 7T 
where Q m and fl^ are the density parameters due to non relativistic matter and cosmological 
constant, respectively, at the time of observation, a and b are 5 the same as A and B with Si 
replaced by Si in their definitions, equations (24) and (12). 

Assuming that the initial density profile has a single, spherically symmetric peak and 
assuming a power law for the initial linear correlation function at the cluster scale, we can 
replace 6(x) by |||, 

£(*) = (^) 3 +", (32) 

where n is the index of the density power spectrum (Peebles 1981) and r 0i is the correlation 
length. This gives 

| = (^ + (3 + O.- "a) (33 ) 

where all the quantities are evaluated at the age of the observed cluster. Note that for models 
of interest the physically plausible range for n is (—3, 0). One can see that the second term 

2 

in equation (33) is indeed proportional to ts and so for an open universe it dominates for 
large time (after curvature domination). Noting that in equation (25), the temperature is 
proportional to £"-f , so, when the second term dominates, the T-M relation will no longer 
evolve with time. This indicates that in an open universe the cluster formation freezes out 
after a certain time. The presence of freeze-out time is independent of the central peak 
approximation since the ratio b/a only depends on the statistics of the initial fluctuations at 
high redshifts where there is very little dependence on cosmology. It is interesting to note 



5 a should not to be confused with the cosmological scale factor. 
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that in the case n = —3, the ratio has no dependence on cosmology and there is no freeze-out 
even in low density universes. This is an interesting case where linear theory does not apply, 
because all scales become nonlinear at the same time and the universe is inhomogeneous on 
all scales. 

Voit (2000) uses a different method to obtain exactly the same result. As we argue next, 
both treatments ignore cluster mergers. 



3.3. B/A: Multiple Peaks and Scatter 

The single peak approximation discussed in §3.2 ignores the presence of other peaks 
in the initial density distribution. In hierarchical structure formation models, the mass of 
a cluster grows with time through mergers as well as accretion. This means that multiple 
peaks may be present within Ri and suppress the effect of the central peak. 

Assuming Gaussian statistics for initial density fluctuations, we can find the statistics 
of b/a. Note that using equation (24), we can fix the value A (and hence a) for a given mass 
and virialization time. So the problem reduces to finding the statistics of b (or B) for fixed 
a. Under the assumption of a power law spectrum (see Appendix A) calculations give 



<b> 4(1 -n) 



a (ra-5)(n-2)' 

with 



(34) 



a 7 16tt2 n/2 r n + 3 ,1,%^ / oc x 

M= (5-n)(2-n) [ n(7-^)(n-3) ]2( ^ ) * ' (35) 

which, inserting into equation (31), yields 

B 4(1 -n) r , 71(71 + 3)., _ _ w Ht .2. . . 

<A >= (n - 5)(n - 2) [1 " " °~ ~ (36) 

^ = r(n)(-^)-^D-Ht)(^)i, (37) 

where D(t) is the growth factor of linear perturbations, normalized to (1 + z)^ 1 for large 
redshift, and 

M 0L = ^Wol, (38) 



- 12 - 



where £c(r) is the linearly evolved correlation function at the present time with r L being 
the correlation length, and 

20 x2-"/ 2 n + 3 i 

r(n)= (5-n)(2-n) [ n(7-^)(^-3) ] ' (39) 

Fig 2 shows the dependence of B/A on n. The upper dotted curve shows the result 
of the single peak approximation (equation 33), while three lower curves show the multiple 
peak calculation described above (equation 36) and its ±lcr dispersion (equation 37). All 
the curves are for an Einstein-de Sitter universe and the dispersion is calculated for mass 
IOMql. Numerically, for r L = 5/i~ 1 Mpc we have Mql = 1.4 x lO u Q m h~ 1 M Q , resulting in 

10M 0L = 4.3 x 10 14 /T 1 M for Vt m = 0.3. We note that n 3 the density distribution is 

dominated by the central peak and corresponds to the top-hat case, and so two methods give 
similar results. Interestingly, as n approaches zero, small peaks dominate and the distribution 
becomes close to homogeneous (top hat) on large scales. This implies that clusters undergo a 
large number of mergers for large values of n. Interestingly, in this case the ratio B/A again 
approaches 2/5, the value for the top-hat case. We will use the multiple-peak approximation 
in our subsequent calculations. 

It is worth mentioning that the ratio of the cosmology dependent term in the average 
value of < B/A > to the constant term, in the multiple-peak calculation, is small. For exam- 
ple, for Q m = 0.3 and n = —1.5, this ratio is about 0.07. This implies that the freeze-out time 
is large comparing to the current age of the universe for feasible open cosmological models 
and consequently < B/A > largely is determined by the spectral index of the underlying 
linear power spectrum n. 



1 



0.8 




-2.5 -2 -1.5 -1 -0.5 

n 

Fig. 2. — The solid curve shows B/A vs. the power index and dashed curves are its ±ler 
variance for M — 10 M 0L (see the text) for an Einstein-de Sitter universe. The dotted curve 
is the result of single peak approximation (equation 33). 
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3.4. v. Surface Term and its Dependence on the final equilibrium density 

profile 

As discussed in §2, corrections to the virial relation due to finite surface pressure changes 
the M-T relation (equation 25). Shapiro, Iliev, & Raga (1999) and Iliev & Shapiro (2001) 
have previously also taken into account the surface pressure term in their treatment of the 
trancated isothermal sphere equilibrium structure with a top-hat initial density perturbation 
and have found results in good agreement with simulations. Voit (2000) uses NFW profile 
for the final density distribution to constrain the extra factor and finds that for typical 
concentration parameters c = r 2 oo/r s ~ 5, [y + — 1) is ~ 2. We will investigate this 

correction, u, for a given concentration parameter c. Let us assume that the density profile 
is given by: 

p(r) = p,f(r/r a ), (40) 

where p s is a characteristic density, r s the scale radius, and / is the density profile. For the 
NFW profile 

W(X) = ^~y 2 , (41) 

Moore et al. (2000) have used simulations with higher resolutions to show that the central 
density profile is steeper than the one already probed by low-resolution simulations such as 
those used by NFW, yielding the Moore profile 

= ^a + s") - (42) 

For a given f(x), the mass of the cluster is: 

M = 4irp s r 3 s g(x), g(x) = f f(x)x 2 dx. (43) 

Jo 

The gravitational energy of the cluster is given by: 

U=-J 9M^L = -lQn 2 Gplr 5 s £ f(x)g(x)xdx. (44) 
To find the surface pressure, we integrate the equation of hydrostatic equilibrium 6 

VP = pg=-^. (45) 



6 This is of course valid in the case of isotropic velocity dispersion profile. As an approximation, we are 
going to neglect any correction due to this possible anisotropy 
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This leads to 

fOO 



f(x)g(x)x- 2 dx, (46) 
which, by the definition of v (equation 15), gives: 

, x ZPextV c 3 f°° f(x)g(x)x- 2 dx , x 

Kc, /) = — = Jc rc y ; . — . (47) 

u Jo f{x)g{x)xdx 

Note that u(c, f) is a function of both c and density profile /. 
If we define Q(c, f) as 

equation (25) can be written as: 

kT = (JpL)(2nGHM)lQ(c, f). (49) 

or in numerical terms: 

A:T = (6.62keV)g(c, / )( iol5/ ^ 1MQ ) 2 / 3 (50) 

for P S pec = 0.9. Note that with this choice of f3 speci the definition of Q is equivalent to that 
of Q in equation (26). 
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Fig (3) shows the average value of y for different cosmologies, using equations 36 and 
51, for n = —1.5. 
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Fig. 4. — Normalization Q vs. y, comparison of different results. The solid and dashed 
curves are for the Moore and NFW profiles respectively. For description of different data 
points see Table 1. The horizontal position of error bars is arbitrary. 
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Since all three parameters, u(c, f), y(c, f) and Q(c, f), are functions of both c and /, one 
can express any one of them as a function of another, for fixed /. Fig (4) shows normalization 
factor Q as functions of y. The dashed curves are for the NFW profile and the solid curves 
for the Moore et al. (2000) profile. The x's with error bars show various simulation and 
observational results (see Table 1). Note that even if we relax the conservation of energy, one 
can still use Q(y) to find the T-M relation, using the value y obtained from its definition, 
equation (49), with the corrected energy. 

An important feature of the behavior of Q(y) is the presence of a minimum in or close 
to the region of physical interest. As a result, Q(y) has very weak dependence on the history 
of the cluster, for example, the largest variation in Q is about 3%. This is probably why 
simulations do not show significant cosmology dependence ~ 5% (e.g. EMN, Mathiesen 



Another way of stating this property is that the heat capacity of the cluster is very 
small. It is well known that the heat capacity of gravitationally bound systems like stars is 
negative. Yet we know that if non-gravitating gas is bound by an external pressure, its heat 
capacity is positive. In the case of clusters, the interplay of external accretion pressure and 
gravitational binding energy causes it to vanish. It is only after the freeze-out time in an 
extreme open universe (low Q m , no cosmological constant) where the heat capacity becomes 
negative similar to an ordinary gravitationally bound system. 



We point out that, by using the conservation of energy, one can constrain the concentra- 
tion parameter and subsequently the surface correction for a given density profile. A typical 
density profile is specified by two parameters: a characteristic density p s and scale radius r s . 
If we know mass (e.g. M 2 oo) and total energy of the cluster, we can fix these two parameters. 
The concentration parameter is then fixed by p s and the critical density of the universe. To 
show the precedure, let us re-derive the T-M relation for a known density profile. Combining 
equations (16) and (19) gives: 



Note that equation (51) only depends on the properties of the virialized cluster and is 
independent of its history. Defining y as 



2000). 



3.5. 



Concentration Parameter c 




(51) 



y=^-(2«GMH )-^, 



(52) 
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equation (51) reduces to: 



Wspec V ~ V f 

Comparing this result with equation (25), we see that 

y = (54) 

only if the energy is conserved [i.e., assuming Ef in equation (49) is equal to Ej\. On the 
other hand, by combining equation (52) with equations (43), (44) and (47) and the virial 
theorem (14), y can be written as a function of c, for a fixed density profile /: 

vie f) = A ' ::i[l ''K^!??' 1 '^''- (55) 



c 



where we have assumed the boundary of the virialized region to be the radius at which the 
average density is A c times the critical density of the universe (which is usually chosen to be 
200), and v is a function of c and / in equation (47). Equation (55) fixes the concentration 
parameter c as a function of y for a fixed density profile /, which in turn is determined by 
equation (54). 

The concentration parameter is fixed by the cosmology (y parameter) as shown in 
Fig. (5). This relation can be well fit by: 

log 10 c= -0.17+ 1.2 y, (56) 

for NFW profile, accurate to 5% in the range < y < 1. 

Let us now consider the evolution of c. We know that in an expanding universe Q m 
decreases with time. Comparing with Fig. (4) we see that in a flat ACDM universe y is 
decreasing with time, while in an open/Einstein-de Sitter universe, it is almost constant. 
Then equation (56) implies that c is a decreasing function of time (increasing function of 
redshift) in a ACDM universe, while it does not significantly evolve in an OCDM universe. 
As an example, the concentration parameter in an Einstein-de Sitter universe is about 40% 
larger than that of a flat ACDM universe with Q m = 0.25. This is consistent with the NFW 
results who find an increase of about 35% for c as a function of mass in units of non-linear 
mass scale. 
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Fig. 6. — Predicted concentration parameter (n = —1.5 vs. mass for Thomas et al. 2001 
simulation of an Einstein-de Sitter cosmology. The shaded area indicates the 68% likelihood 
region. The dashed line is the best fit to their result. 
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NFW simulations show a weak dependence of the concentration parameter on mass 
c oc M~ 0A . We see that our concentration parameter does not depend on mass. However its 
scatter is larger for small masses and so is marginally consistent with the simulation results 
(see Fig 6). Assuming that this discrepancy is only a consequence of non-spherical shape of 
original proto-cluster, in the next section, we attempt to modify the value of y to match the 
simulation results. 



3.6. Corrections for Initial Non-Sphericity 

In this section we try to incorporate the effects due to the non-spherical shape of the 
initial proto-cluster into our formalism. Unlike previous sections, the calculations of this 
section are not very rigorous and should be considered as an estimate of the actual correc- 
tions. In particular, these approximations lose accuracy if there are large deviations from 
sphericity which, as we see, is the case for low mass end of the M-T diagram. 

We are going to assume that non-sphericity comes in through a modifying factor 1 + M 
that only depends on the initial geometry of the collapsing domain, 

y N = y{l+N) ) (57) 

where y^- is the modified value of y. Next, let us assume that Ri(9,ip) is the distance of 
the surface of our collapsing domain from its center. We can expand its deviation from the 
average in terms of spherical harmonics Yi m (9,<p), 

5R i (9,<p) = J2 a i™Yim(0,<p). (58) 

l,m 

If we try to write down a perturbative expansion for A/", the lowest order terms will be 
quadratic, since there is no rotationally invariant first order term. Moreover, having in mind 
that the gravitational dynamics is dominated by the large scale structure of the object, as 
an approximation, we are going to keep the lowest I value. Since I — 1 is only a translation 
of the sphere and does not change its geometry, the lowest non-vanishing multipoles are for 
I = 2, and the only rotation invariant expression is 

2 

m=-2 

where we absorbed any constant factor in the definition of c^m's. The next simplifying as- 
sumption is that a 2m 's are Gaussian variables, with amplitudes proportional to the amplitude 
of the density fluctuations at the cluster mass scale. Mathematically, this is motivated by 
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the fact that the concentration parameter predicted by simulations is closer to our prediction 
for spherical proto-clusters, at large mass end, where amplitude is smaller. The physical mo- 
tivation is that since the density fluctuations decrease with scale, more massive clusters tend 
to deviate less from sphericity Choosing Gaussian statistics for a 2 mS is only a simplifying 
assumption to carry out the calculations. Then it is easy to see that 

AAf 2 = <Af > 2 , (60) 
25 

and then using the definition of Af, assuming that it is a small correction we get 

(^) 2 --(l-^) 2 + A 2 . (61) 

Note from equation (54), Ay = (Ht)~ 2 ^(AB / A) and AB/A is given in equation (37). We 
have also assumed that Af and y are statistically independent variables. In the next step, 
we define the amplitude of Af 

AB 

< Af >^ u(—) 2 , uj^Q4 (62) 

The numerical value of uo is fixed by plugging y^ into equation (56) to get the modified 
concentration parameter and comparing this result with the simulations of Thomas et al. 
(Fig 6). Fig 7 shows the modified concentration parameter as a function of mass in an 
Einstein-de Sitter universe. We see that the introduction of non-sphericity results in the 
cluster concentration parameter being a decreasing function of cluster mass with a scatter 
that also decreases with mass, in accord with simulations. Fig 8 shows the same comparison 
for a ACDM cosmology with Eke et al. (2001) fitting formula. We see that, although u was 
obtained by fitting the Einstein-de Sitter simulations, our prediction is marginally consistent 
with the ACDM simulations as well. 7 



7 Ekc ct al. (2001) fitting formula was made for a lower mass range and its systematic error at the cluster 
mass range is indeed comparable to its difference with our prediction. 
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Fig. 7. — Similar to Fig 6 but with non-spherical correction. The solid curve is fitted to 
match the simulation result. 



-26- 




13 13.5 14 14.5 15 

log(M 200 h/M o ) 

Fig. 8. — The modified concentration parameter vs. mass for a ACDM cosmology with 
Q m = 0.3, Q\ = 0.7, n = —1.5 and a 8 = 0.9. The dashed line is the Eke et al. (2001) fitting 
formula with C a = 23 and the shaded area is the 68% likelihood region. 
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3.7. Scatter in M-T relation 

Since y is linear in B/A, it also has a Gaussian probability distribution function (PDF): 

P{y)dy = dy exp[- (y ~^ 2 ], (63) 

where y and Ay are related to equations (36) and (37), by equation (54). As we mentioned 
above, the variation of Q for the average values of y is negligible. However, the scatter in the 
value of y can be large, especially in the low mass end (see equation 37) and so the scatter 
in Q might become significant. In what follows we only consider the NFW profile, since it 
is extensively considered in the literature. We find that the behavior of Q can be fitted by: 

Q( y y = Ql + N(y-y ) 2 , (64) 

where 

Qo = 1.11, AT = 1.8, y = 0.538. (65) 

The error in this fitting formula is less than 3% in the range — 1 < y < 2. Inserting this into 
equation (63) leads to the PDF of Q: 

P(Q)dQ = , 2QdQ -. 



exp[-^±g^] coshpg# V^l- (66) 

Fig 9 shows three different examples of the PDF obtained here. It is clear that the 
scatter in Q is asymmetric. In fact, since the average value of y is close to the minimum of 
equation (57), the scatter in y shifts the average value of Q upwards systematically. In the 
limit of large Ay, where this shift is significant, P(Q) is approximately: 

PmdQ ^ exph w 1 ' 9(0 - w ' (67) 



where 6 is the Heaviside step function. 





1 1.2 1.4 1.6 

Q 

Fig. 9. — Probability distribution function of normalization for three different y and Ay's. 
The PDF vanishes for Q < Q — 1.11. The large value of y = 1.0 is possible if one includes 
non-spherical corrections (§3.6). 
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Although the assumption of Gaussianity is not strictly valid for y^f which includes non- 
spherical corrections, we still can, as an approximation, use the above expressions for the 
PDF by replacing y and Ay by y^ and Ay^. It is also easy to find the average of Q 2 

<Q 2 >=Ql + N[(y-y ) 2 + Ay 2 ]. 

As an example, for y ~ y , while Ay = 0.5 gives only about 16% systematic increase in Q, 
Ay = 1.0 leads to~ 60% increase. 

3.8. Observed and Average Temperatures 

The temperature found in §2.4 is the density weighted temperature of the cluster aver- 
aged over the entire cluster. However, the observed temperature, Tf, can be considered as 
a flux-weighted spectral temperature averaged over a smaller, central region of the cluster. 
The two temperatures may be different, due to presence of inhomogeneities in temperature. 
We use the simulation results by Mathiesen & Evrard (2001) to relate these two tempera- 
tures and refer the reader to their paper for their exact definitions and the discussion of the 
effects which lead to this difference: 

T — Tf[l + (0.22 ± 0.05) log 10 7>( keV) - (0.11 ± 0.03)]. (68) 

This correction changes the mass-temperature relation from M oc T L5 to M oc T for 
the observed X-ray temperatures. We use this correction in converting the observed X-ray 
temperature to virial temperature in Figures 1 and 10-12. 

4. Predictions vs. Observations 
4.1. Power Index 

It is clear from equation (49) that we arrive at the usual M oc T 15 relation that is 
expected from simple scaling arguments and is consistent with the numerical simulations 
(e.g. EMN or Bryan & Norman 1998). On the other hand, the observational /9-model 
mass estimates lead to a steeper power index in the range 1.7 — 1.8. Although originally 
interpreted as an artifact of the /3-model (HMS), the same behavior was seen for masses 
estimated from resolved temperature profile (FRB). FRB carefully analyzed the data and 
interpret this behavior as a bent in the M-T at low temperatures. This was confirmed by 
Xu, Jin & Wu (2001) who found the break at T x = 3 - 4 keV. 
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As discussed in §3.7, the asymmetric scatter in Q introduces a systematic shift in the 
M-T relation. For large values of Ay all of the temperatures are larger than the value given 
by the scaling relation (54) for average value of y. This scatter increases for smaller masses 
(see equation 37), hence smaller clusters are hotter than the scaling prediction. As a result, 
the M-T relation becomes steeper in the low mass range as Q increases, while the intrinsic 
scatter of the data is also getting larger. Indeed, increased scatter is also observed in the 
FRB data(Fig 1 and Fig's 10-12 to compare with our prediction) but they interpret it as the 
effect of different formation redshifts. We will address this interpretation in §5. 



4.2. Normalization 

As we discussed in §3.4, our normalization (i.e., Q in equation 50) is rather stable with 
respect to variations of cosmology and the equilibrium density profile. Table 1 compares this 
value with various observational and simulation results. 
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Q 


Method 


Reference 


1.12 ± 0.02 

1.21 ±0.06 
1.05 ±0.13 

1.22 ±0.11 
1.32 ±0.17 
1.70 ±0.20 
1.70 ±0.18 


Analytic 
Hydro-Simulation 
Hydro-Simulation 
Hydro-Simulation 
Optical mass estimate 
Resolved temperature profile 
Resolved temperature profile 


This paper 
Thomas et al. 2001 
Bryan & Norman 1998 (BN) 
Evrard et al. 1996 (EMN) 
Horner et al. 1999 (HMS1) 
Horner et al. 1999 (HMS2) 
Finoguenov et al. 2001 (FRB) 



Table 1. Comparison of normalizations from different methods. The last two rows only 

include clusters hotter than 3 keV. 
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We see that our analytical method is consistent with the hydro-simulation results, in- 
dicating validity of our method, since both have a similar physics input and the value of 
fispec used here was obtained from simulations. On the other hand, X-ray mass estimates 
lead to normalizations about 50% higher than our result and simulations. The result of 
the optical mass estimates quoted is marginally consistent with our result. Assuming that 
this is true would imply that optical masses are systematically higher than X-ray masses 
by ~ 80%. Aaron et al. (1999) have compared optical and X-ray masses for a sample of 
14 clusters and found, on the contrary, a systematic difference less than 10%. The optical 
masses used in the work of Aaron et al. (1999) were all derived by CNOC group (Carlberg 
et al. 1996) while the Girardi et al. 1998 data, used in HMS analysis quoted above, was 
compiled from different sources and has larger scatter and unknown systematic errors. In 
fact HMS excluded a number of outliers to get the correct slope and original data had even 
larger scatter (see their Fig 1). Therefore it may be that the systematic error in the optical 
result of Table 1, be much larger and so in agreement with other observations. 

As discussed in §3.1, one possible source for difference between theoretical and obser- 
vational normalizations is that the values for f3 spec are different in the two cases due to 
systematic selection effects. Also, intriguingly, Bryan & Norman (1998) show that there is a 
systematic increase in the obtained value of f3 spec by increasing the resolution of the simula- 
tions. Whether this is a significant and/or real effect for even higher resolutions is not clear 
to us. 

However, the fact that the slope is unchanged indicates that the missing process is, 
probably, happening at small scales and so relates the intermediate-scale temperature to 
the small-scale flux-weighted spectral temperature by a constant factor, independent of the 
large-scale structure of the cluster. In this case the actual value of (3 spec must be ~ 0.6. 

Figures 10-12 show the prediction of our model, shifted downwards to fit the observa- 
tional data in the massive end, versus the observational data of FRB using resolved temper- 
ature profile and corrected as discussed in §3.8. The correction due to initial non-sphericity 
(§3.6) is included in the theoretical plot. The value of a s which enters AB/A (equation 37) 
through M 0L , is fixed by cluster abundance observations (e.g. Bahcal & Fan 1998). 
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Fig. 10. — Our result with shifted normalization for an Einstein-de Sitter universe vs. FRB 
data. The shaded area indicates the 68% confidence level region (§3.7). 



-34- 




0.5 1 

log(T(keV)) 

Fig. 11. — Our result with shifted normalization for an OCDM universe vs. FRB data. The 
shaded area indicates the 68% confidence level region (§3.7). 
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We see that while an Einstein-de Sitter cosmology under-estimates the scatter in the 
low mass end, a typical low density OCDM cosmology overestimates it. On the other hand, 
a typical ACDM cosmology is consistent with the observed scatter. Interestingly, this is 
consistent with various other methods, in particular CMB+SNe la result which point to a 
low-density flat universe (de Bernardis et al. 2000; Balbi et al. 2000; Riess et al. 1998). 



4.3. Evolution of M-T Relation 

As we discussed in §3.4, the value of our normalization has a weak dependence on 
cosmology. Going back to equation (49), we see that the time dependence of the M-T 
relation is simply: M oc /7~ 1 T 1 ' 5 . Assuming a constant value of /9 spec , this formula can be 
potentially used to measure the value of H at high redshifts and so constrain the cosmology. 

Schindler (1999) has compiled a sample of 11 high redshift clusters (0.3 < z < 1.1) 
from the literature with measured isothermal /9-model masses. In these estimates, the gas is 
assumed to be isothermal and have the density profile: 

P,W=P 3 (0)(l + (f) 2 )- 3 ^/ 2 (69) 

' c 

Then, the mass in overdensity A c is given by: 

m - 0" /2 < t^' 3/2 < 7o > 

Comparing this with equation (49), and neglecting the difference between virial and X-ray 
temperatures, we get: 

Q ~ (37r)- 2 / 3 (^ ec //3 /it )Ay 3 (71) 

In the last two equations, we have ignored r c with respect to the radius of the virialized 
region, which introduces less than 3% error. This allows us to find the normalization Q from 
the value of (3f it (independent of cosmology in this case). 

Assuming f3 spec = 0.9, equation (71) gives the value of Q for a given f3f it . Fig 13 shows 
the value of Q versus redshift for Schindler (1999) sample and also FRB resolved temperature 
method for low redshift clusters. This result is consistent with no redshift dependence and 
the best fit is: 



logQ = 0.23 ± OM(systematic) ± 0M(random) + (0.09 ± 0.04)z. 



(72) 
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Fig. 13. — Evolution of the normalization Q with redshift. The data points are from Schindler 
(1999) high-z sample and FRB low-z normalization at z — 0. The line is the best fit. 
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The combination of this result and equation (50) gives 



kT= (11.2 ± 1.1 keV)e( a21±0 09 ^( 



M 




(73) 



© 



which relates X-ray temperature of galaxy clusters to their masses that does not depend 
on the theoretical uncertainties with regard to the normalization coefficient of the M-T 
relation, in this range of redshifts. The systematic error in this result is less than 5% while 
the random scatter can be as large as 20%. This result is valid for T x > 4keV since below 
this temperature the systematic shift due to random scatter becomes significant (§3.7). It is 
easy to see that this threshold moves to lower temperatures in high redshifts in our formalism. 

Note that the more realistic interpretation of a possible evolution in observed Q, ob- 
tained above, is that in equation (71), Q remains constant (Fig 5) and, instead, f3 spec varies 
with time. However there will be virtually no difference with respect to the M-T relation 
and, moreover, the weak redshift dependence in equation (73) shows that f3 spec is indeed 
almost constant. 



In this section we discuss the validity of different approximations which were adopted 
throughout this paper. 

In the calculation of initial and final energy of the cluster, we ignored any contribution 
from the vacuum energy. In fact we know that cosmological constant in the Newtonian limit 
can be considered as ordinary matter with a constant density. If the cosmological constant 
does not change with time, then its effect can be considered as a conservative force and so 
energy is conserved. However, in both initial state and final equilibrium state of the cluster 
the density of cosmological constant is much smaller than the density of matter and hence its 
contribution is negligible. This does not hold for quintessential models of the vacuum energy 
since A changes with time and energy is not conserved. This may be used as a potential 
method to distinguish these models from a simple cosmological constant. 

Let us make a simple estimate of the importance of this effect. We expect the relative 
contribution of a varying cosmological constant be maximum when the density of the proto- 
cluster is minimum. This happens at the turn-over radius which is almost twice r v , the 
virial radius, in the top-hat approximation. Then the contribution to the energy due to the 
vacuum energy would be: 



5. 



Discussion 




5E 



E 



200 p c (t ) 



8 v A(t /2) 




8 /(W2) 
200 A H 2 (t ) 



)Ov(V2) ~ 0.1. (74) 
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This gives about 10% correction to y or about 15% correction to c. We see that still the 
effect is small but, in principle observable if we have a large sample of clusters with measured 
concentration parameters. 

A systematic error in our results might have been introduced by replacing the density 
distribution by an averaged radial profile which under-estimates the magnitude of gravita- 
tional energy and so the temperature. Also, as shown by Thomas et al. (2000), most of the 
clusters in simulations are either steeper or shallower than the NFW profile at large radii. 
Since our normalization is consistent with simulation results, we think that these effects do 
not significantly alter our predictions. 

Let us now compare our results with that of Voit (2000), who made the first such 
analytic calculation (but we were not aware of his elegant work until the current work was 
near completion) which is not based on the top-hat initial density perturbation and uses the 
same ingredients to obtain M-T relation. First of all, as noted in §3.3, his result, which is 
equivalent to the central peak approximation to find B/A, ignores the possibility of mergers. 
As we see in Fig (2) the value of B/A and so y is about 50% larger in the single peak case 
than the multiple peak case. Although it does not change the normalization very much, it 
overestimates the concentration parameter (Fig (5)) in our formalism. However, as pointed 
out by the referee, our formaism for finding the value of c is based on the assumption of 
isotropic velocity dispersion (equation (45)) and hence may not be directly compared to Voit 
(2000) who does not make such assumption. Also Voit (2000) has neglected the cosmology 
dependence of the surface correction, which gives a false cosmology dependence to the 
normalization. This is inconsistent with hydro-simulations (e.g. EMN, Mathiesen 2000) and 
does not give the systematic shift in the lower mass end, if we assume it does not have a 
non-gravitational origin. 

Finally, we comment on the interpretation of the scatter/bent in the lower mass end 
as being merely due to different formation redshifts which is suggested by FRB. Although 
different formation redshifts can certainly produce this scatter/bent, it is not possible to 
distinguish it from scatter in initial energy of the cluster or its initial non-sphericity, in our 
formalism. On the other hand, the FRB prescription, which assumes constant temperature 
after the formation time, is not strictly true because of the on-going accretion of matter even 
after the cluster is formed. So, as argued by Mathiesen 2000 using simulation results, the 
formation time might not be an important factor, whereas the effect can be produced by 
scatter in initial conditions of the proto-cluster. 
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6. Conclusion 

We combine conservation of energy with the virial theorem to derive the mass-temperature 
relation of the clusters of galaxies and obtain the following results: 

• Simple spherical model gives the usual relation T = CM 2 / 3 , with the normalization 
factor C being in excellent agreement with hydro simulations. However, both our 
normalization and that from hydro simulations are about 50% higher than the X-ray 
mass estimate results. This is probably due to our poor understanding of the history 
of the cluster gas. 

• Non-sphericity introduces an asymmetric, mass-dependent scatter (the lower the mass, 
the larger the scatter) for the M — T relation thus alters the slope at the low mass 
end (T ~ 3keV). We can reproduce the recently observed scatter/bent in the M-T 
relation in the lower mass end for a low density ACDM cosmology, while Einstein- de 
Sitter/OCDM cosmologies under/overestimate this scatter/bent. We conclude that the 
behavior at the low mass end of the M-T digram can be used to constrain cosmological 
models. 

• We point out that the concentration parameter of the cluster and its scatter can be 
determined by our formalism. The concentration parameter determined using this 
method is marginally consistent with simulation results, which provides a way to find 
non-spherical corrections to M-T relation by fitting our concentration parameter to the 
simulation results. 

• Our normalization has a very weak dependence on cosmology and formation history. 
This is consistent with simulation results. 

• We find mass-temperature relation (73), for clusters of galaxies, based on the observa- 
tions calibrated by our formalism, which can be used to find masses of galaxy clusters 
from their X-ray temperature in the range of redshift < z < 1.1 with the accuracy of 
20%. This is a powerful tool to find the evolution of mass function of clusters, using 
their temperature function. 

This research is supported in part by grants NAG5-8365. N.A. wishes to thank Ian 
dell' Antonio, Licia Verde and Eiichiro Komatsu for useful discussions. 
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Assuming a Gaussian statistics for the linear density field, The Probability Distribution 
Function for a and b, takes a Gaussian form: 



A. Statistics of b/a 



a and b are defined as: 




(Al) 



P(a, b) da db 




(A2) 



This preprint was prepared with the AAS IATeX macros v5.0. 
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a < a 2 > 



Then, for fixed a, we have: 

<6> <a6> (A3) 

*» = V^> < A4 > 

To find the quadratic moments of a and b, we assume a power-law linear correlation function: 

fc(r) =< ^(x)^(x + r) >= (^) 3+n . (A5) 

The moments become: 

<a 2 > = 87T 2 (^) 3+ "F o(n), (A6) 

Hi 

<a b> = 8 7r 2 (^) 3 +"(F 00 (n)-F 02 (n)), (A7) 
<b * > = 8 7r 2 (^) 3+ "(F 00 (n)-2F 02 (n) + F 22 (n)), (A8) 

where F m i(n) is defined as: 

F ml (n) = -L / <4| Xl - x.r^^X!^, (A9) 

8^ J 

where the integral is taken inside the unit sphere. Taking the angular parts of the integral, 
this reduces to: 

F ml = — *— / / dx 1 dx 2 x? +1 x l 2 +1 (\x 1 - x 2 \~ 1 - n - \ Xl + x 2 \- l ~ n ). (A10) 
n + 1 Jo Jo 

Then, taking this integral for the relevant values of m and I, and inserting the result into 
(A6-A8) and subsequently (A3-A4) gives: 

±*±=, 4(1 ~ n) „ (All) 

a (n-5)(n-2)' V ; 

with 

A6 = 167r2 ~" /2 [ ri±3 (A12) 

(5-n)(2-n) [ n(7-n)(n-3) J 1 ; 



